!.....7..0.........0.........0.........0.........0.........0.........012
      PROGRAM RSLICES           
!.....7..0.........0.........0.........0.........0.........0.........012
!
#include "param.nml"
      dimension th_mlt(l)
      call analytic(th_mlt)
      stop 'main'
      end   
!
!
      subroutine the_prof_mlt(th_mlt)
!
!  Stratification using mixinglength relation for the
!  convective Flux
!
!  dT/dx=g(x)/(cv*(gamma-1))*Grad
!  dlnrho/dx=g(x)/(cv*(gamma-1))*(1-Grad)/T
!
!  Grad_cool = 0
!  Grad_conv = Grad_ad + k*(Fmlt/rho*c**3)**2/3
!  Grad_rad = 1/(mpoly_ml0+1)
!
!  11-March-11/Gustavo: Based on Boris routine and
!  (Brandenburg et al. AN, 326, 681, 2005)

#include "param.nml"

      dimension pm0(l),tm0(l),z(l),rad(l),tt(l),rho(l)
      dimension rlnrhom(l),tempm(l),th_mlt(l)
      real rmin,rmax,dr
      real zm,zbot,ztop, Fmlt,rlnrho,ss,rlnTT,sluminosity
      real rbot, rt_old, rt_new, rb_old, rb_new,gamma,gamma_m1,cap
      real rhotop, rhobot,crit
      real rho_ic, cs2top_ic,rho0,rlnrho0,tt0,lntt0,th,th0,thi
      real z1, z2, width_ic,solar_radii
      real temp00,rho00,th00
      integer ipz,k,iter,iz
!
!      solar_radii=6.96e8
      solar_radii=1.
      rmin=0.6*solar_radii
      rmax=0.96*solar_radii
!
      dr=(rmax-rmin)/(l-1)
      do 11 k=1,l
   11 rad(k)=rmin+dr*(k-1)

      ipz=0
      gamma=5./3.
      gamma_m1=gamma-1.
      cap=1.-1/gamma
! normalization factors
      temp00=3.4976e6
      rho00=370.33
      th00=temp00

!      z1=0.71*solar_radii
!      z2=1.*solar_radii
!      width_ic=0.005*solar_radii
!      rho0_ic = 190.9
!      cs2top_ic = 213422.05*gamma_m1
!      sluminosity=3.846e11
      z1=0.71
      z2=1.2
      width_ic=0.02
      rho0_ic = 10.
      cs2top_ic = 0.078432
      sluminosity=0.002

      ztop=rad(l)
      zbot=z1
      rbot = rho0_ic
      rt_old=.1*rbot
      rt_new=.12*rbot
!
      Fmlt = sluminosity/(4.*3.1416*rad(1)**2)
!
!  Need to iterate for rhobot=1.
!  Produce first estimate.
      rhotop = rt_old
      cs2top = cs2top_ic
      call strat_MLT(rhotop,Fmlt,rlnrhom,tempm,rhobot,z1,z2,rad)
      rb_old=rhobot
!
! Next estimate
!
      rhotop=rt_new
      call strat_MLT(rhotop,Fmlt,rlnrhom,tempm,rhobot,z1,z2,rad)
      rb_new=rhobot
!
      do 10 iter=1,10
!
!  New estimate.
!
         rhotop=rt_old+(rt_new-rt_old)/(rb_new-rb_old)*(rbot-rb_old)
!
!  Check convergence.
!
         crit=abs(rhotop/rt_new-1.)
         if (crit<=1e-6) goto 20
!
         call strat_MLT(rhotop,Fmlt,rlnrhom,tempm,rhobot,z1,z2,rad)
!
!  Update new estimates.
!
        rt_old=rt_new
        rb_old=rb_new
        rt_new=rhotop
        rb_new=rhobot
10    continue
20    print*,'completed: rhotop,crit=',rhotop,crit,iter
!
      rlnrho0 = rlnrhom(l)
      rho0 = exp(rlnrho0)
      tt0 = tempm(l)
!      
      open(1,file="strat.dat")
!
      do k=1,l
        iz=k+ipz*l
        zm=rad(1)+(iz-1)*dr
        rlnrho=rlnrhom(l-iz+1)
        rlnTT=log(tempm(l-iz+1))
        rho(k)=rho00*exp(rlnrho)/rho0
        tt(k)=temp00*tempm(l-iz+1)/tt0
      enddo
!
      do k=1,l 
        th=tt(k)*(rho(l-1)*tt(l-1)/(tt(k)*rho(k)))**cap
        th_mlt(k)=th
        write(1,'(F5.4,x,F13.3,x,F8.4,x,F13.3)'),rad(k),&
                  tt(k),rho(k), th_mlt(k)
        print*,rad(k), tempm(k), rho(k), th_mlt(k)
      enddo
!
      close(1)
!
      return
      end
!
      subroutine strat_MLT(rhotop,solar_flux,rlnrhom,tempm,&
                           rhobot,z1,z2,rad)
!
#include "param.nml"
!
      dimension pm0(l),tm0(l),z(l),rad(l)
      dimension rlnrhom(l),tempm(l),zz(l)
      real rhotop, zm, ztop, dlnrho, dtemp,k_mlt,dr
      real solar_flux, rlnrhobot, rhobot, z1,z2
      real del, delad, fr_frac, fc_frac, fc, polyad
      real cs2top,mpoly0,mpoly1,mpoly2,gamma,gamma_m1
      real rlnr_tmp, etmp,solar_mass, GG, gravz,solar_radii
      integer iz,izm1,k,nbot1,nbot2
!
      dr=rad(2)-rad(1)
      k_mlt = 2.5
      gamma=5./3.
      gamma_m1=gamma-1.
      mpoly0=2.
      mpoly1=1.4
      mpoly2=mpoly1
!      GG=6.674e-11   !N(m/kg)^2
!      solar_mass=1.988e30     !kg
!      cs2top = 213422.05*gamma_m1
!      solar_radii=6.96e8
      GG=3.
      solar_mass=1.
      cs2top = 0.078432
      solar_radii=1.
!
!  Initial values at the top.
!
      rlnrhom(1)=alog(rhotop)
      tempm(1)=cs2top/gamma_m1
!
      ztop=rad(l)
      zz(1)=ztop
!
      polyad=1./gamma_m1
      delad=1.-1./gamma
      fr_frac=delad*(mpoly1+1.)
      fc_frac=1.-fr_frac
      fc=fc_frac*solar_flux
!
      do iz=2,l
        izm1=iz-1
        zm=ztop-(iz-1)*dr
        zz(iz)=zm
        gravz=-GG*solar_mass/(zm**2)
        if (zm<=z1) then
!
!  radiative zone=polytropic stratification
!
          del=1./(mpoly0+1.)
        else
          if (zm<z2) then
!
!  convective zone=mixing-length stratification
!
            rlnr_tmp=rlnrhom(iz-1)
            etmp=exp(rlnr_tmp)*(gamma_m1*tempm(iz-1))**1.5
            del=delad+k_mlt*(fc/etmp)**.6666667
!
          else
!
!  upper zone=isothermal stratification
!
            del=0.
          endif
        endif
        dtemp=gamma*polyad*gravz*del
        dlnrho=gamma*polyad*gravz*(1.-del)/tempm(iz-1)
        tempm(iz)=tempm(iz-1)-dtemp*dr
        rlnrhom(iz)=rlnrhom(iz-1)-dlnrho*dr
      enddo
!
!  Find the value of rhobot.
!
      do iz=1,l
        if (zz(iz)<=z1) exit
      enddo
      nbot1=iz-1
      nbot2=iz
!
!  Interpolate.
!
      rlnrhobot=rlnrhom(nbot1)+(rlnrhom(nbot2)-rlnrhom(nbot1))/&
               (zz(nbot2)-zz(nbot1))*(z1-zz(nbot1))
      rhobot=exp(rlnrhobot)
!
      return
      end

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
      subroutine tripoly(th_mlt)
!
#include "param.nml"
!
      dimension pm0(l),tm0(l),z(l),rad(l),th_mlt(l)
      dimension rlnrhom(l),tempm(l),zz(l),dTdr(l),dTedr(l),rhom(l)
      dimension dlnTdr(l), dlnrdr(l),cs2(l),gravz(l)
      real solar_flux, rlnrhobot, rhobot, z1,z2, width_ic
      real mpoly, mpoly0,mpoly1,mpoly2,gamma,gamma_m1
      real solar_mass, GG, solar_radii,r,th
!
      solar_radii=6.96e8
      rmin=0.60*solar_radii
      rmax=0.96*solar_radii
      z1=0.718*solar_radii
      z2=1.99*solar_radii
      width_ic=0.015*solar_radii
      dr=(rmax-rmin)/(l-1)
!
      ipz=0
      gamma=5./3.
      gamma_m1=gamma-1.
      cap=1.-1/gamma
      Rg=13732.
      cp=2.5*rg
      cv=cp/gamma 
!
!! starting at r=0.62R 
!      temp00=2.97382e+06
!      rho00=420.92
!      th00=2.0868e6
! starting at r=0.60R
      rho00= 183.79604    
      temp00=2129425.9
      th00=3.75236e6
      lbcz=31
!
      mpoly0=3.
      mpoly1=1.49999
      mpoly2=mpoly1

      GG=6.674e-11   !N(m/kg)^2
      solar_mass=1.988e30     !kg
      tempm(lbcz)=temp00
      rlnrhom(lbcz)=alog(rho00)
      rt_new=4.
!      
      do k=1,l
        rad(k)=rmin+dr*(k-1)
        gravz(k)=-GG*solar_mass/(rad(k)**2)
        mpoly=mpoly0-(mpoly0-mpoly1)*0.5*(1.+erf((rad(k)-z1)/width_ic)) &
             -(mpoly1-mpoly2)*0.5*(1.+erf((rad(k)-z2)/width_ic))
        dTdr(k)=gravz(k)/(cv*gamma_m1*(mpoly+1.))
      enddo
!
      do k=lbcz-1,1,-1
        tempm(k)=tempm(k+1)-dTdr(k)*dr
      end do
      do k=lbcz+1,l
        tempm(k)=tempm(k-1)+dTdr(k-1)*dr
      enddo

      do k=1,l
        r=rmin+dr*(k-1)
        dTedr(k)=deriv3(r,rad,tempm,l,1)
        dlnTdr(k)=dTdr(k)/tempm(k)
        dlnrdr(k)=gravz(k)/(tempm(k)*cv*(gamma-1.))-dlnTdr(k)
      enddo
      do k=lbcz-1,1,-1
        rlnrhom(k)=rlnrhom(k+1)-dlnrdr(k+1)*dr
        rhom(k)=exp(rlnrhom(k))
      enddo
     do k=lbcz+1,l
        rlnrhom(k)=rlnrhom(k-1)+dlnrdr(k-1)*dr
        rhom(k)=exp(rlnrhom(k))
     enddo

       rhom(lbcz)=rho00

      open(1,file="strat.dat")
      do k=1,l 
        th=tempm(k)*(rhom(1)*tempm(1)/(tempm(k)*rhom(k)))**cap
        th_mlt(k)=th
        write(1,'(F5.4,x,F13.3,x,F8.4,x,F13.3)'),rad(k)/solar_radii,&
                  tempm(k),rhom(k), th_mlt(k)
        print*,k,rad(k)/solar_radii, tempm(k), rhom(k), th_mlt(k)
      enddo
!
      close(1)
!
      return 
      end

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
      subroutine tripoly2(th_mlt)
!
#include "param.nml"
!
      dimension pm0(l),tm0(l),z(l),rad(l),th_mlt(l)
      dimension rlnrhom(l),tempm(l),zz(l),dTdr(l),dTedr(l),rhom(l)
      dimension dlnTdr(l), dlnrdr(l),cs2(l),gravz(l)
      real solar_flux, rlnrhobot, rhobot, z1,z2, width_ic
      real mpoly, mpoly0,mpoly1,mpoly2,gamma,gamma_m1
      real solar_mass, GG, solar_radii,r,th
!
      solar_radii=6.96e8
      rmin=0.60*solar_radii
      rmax=0.96*solar_radii
      z1=0.718*solar_radii
      z2=1.99*solar_radii
      width_ic=0.005*solar_radii
      dr=(rmax-rmin)/(l-1)
!
      ipz=0
      gamma=5./3.
      gamma_m1=gamma-1.
      cap=1.-1/gamma
      Rg=13732.
      cp=2.5*rg
      cv=cp/gamma 
!
!! starting at r=0.62R 
!      temp00=2.97382e+06
!      rho00=420.92
!      th00=2.0868e6
! starting at r=0.60R
      rho00= 183.79604    
      temp00=2129425.9
      th00= 521.     
!
      mpoly0=3.
      mpoly1=1.499
      mpoly2=mpoly1

      GG=6.674e-11   !N(m/kg)^2
      solar_mass=1.988e30     !kg
      tempm(1)=temp00
      rlnrhom(1)=alog(rho00)
!      
      do k=1,l
        rad(k)=rmin+dr*(k-1)
        gravz(k)=-GG*solar_mass/(rad(k)**2)
        mpoly=mpoly0-(mpoly0-mpoly1)*0.5*(1.+erf((rad(k)-z1)/width_ic)) &
             -(mpoly1-mpoly2)*0.5*(1.+erf((rad(k)-z2)/width_ic))
        dTdr(k)=gravz(k)/(cv*gamma_m1*(mpoly+1.))
      enddo
!
      do k=2,l
        tempm(k)=tempm(k-1)+dTdr(k-1)*dr
      enddo

      do k=1,l
        cs2(k)=cv*tempm(k)*(gamma-1.)*gamma
        r=rmin+dr*(k-1)
        dTedr(k)=deriv3(r,rad,tempm,l,1)
        dlnTdr(k)=dTedr(k)/tempm(k)
        dlnrdr(k)=gravz(k)/(tempm(k)*cv*(gamma-1.))-dlnTdr(k)
      enddo
      do k=2,l
        rlnrhom(k)=rlnrhom(k-1)+dlnrdr(k-1)*dr
        rhom(k)=exp(rlnrhom(k))
      enddo
        rhom(1)=exp(rlnrhom(1))

      open(1,file="strat.dat")
      do k=1,l 
        th=tempm(k)*(rhom(1)*tempm(1)/(tempm(k)*rhom(k)))**cap
        th_mlt(k)=th00*th/temp00
        write(1,'(F5.4,x,F13.3,x,F8.4,x,F13.3)'),rad(k)/solar_radii,&
                  tempm(k),rhom(k), th_mlt(k)
        print*,rad(k)/solar_radii, tempm(k), rhom(k), th_mlt(k)
      enddo
!
      close(1)
!
      return 
      end

      subroutine analytic(th_mlt)
!
#include "param.nml"
!
      dimension pm0(l),tm0(l),z(l),rad(l),th_mlt(l)
      dimension rlnrhom(l),tempm(l),zz(l),dTdr(l),dTedr(l),rhom(l)
      dimension dlnTdr(l), dlnrdr(l),cs2(l),gravz(l),zcr(l)
      real solar_flux, rlnrhobot, rhobot, z1,z2, width_ic
      real mpoly, mpoly0,mpoly1,mpoly2,gamma,gamma_m1
      real solar_mass, GG, solar_radii,r,th,hb,cap,g0
!
      solar_radii=6.96e8
      rmin=0.61*solar_radii
      rmax=0.96*solar_radii
      z1=0.7*solar_radii
      z2=1.96*solar_radii
      width_ic=0.015*solar_radii
      dr=(rmax-rmin)/(l-1)
!
      ipz=0
      gamma=5./3.
      gamma_m1=gamma-1.
      cap=1.-1/gamma
      Rg=13732.
      cp=2.5*rg
      cv=cp/gamma 
      g0=734.428
!
! starting at r=0.60R
      temp00=3.4976e6
      rho00=370.33
      th00=temp00
      
!
      mpoly0=1.8
      mpoly1=1.499995 !1.4999989
      mpoly2=1.499995

      GG=6.674e-11   !N(m/kg)^2
      solar_mass=1.988e30     !kg
      tempm(1)=temp00
      rhom(1)=rho00

      do k=1,l
        rad(k)=rmin+dr*(k-1)
        zcr(k)=dr*k
        mpoly=mpoly0-(mpoly0-mpoly1)*0.5*(1.+erf((rad(k)-z1)/width_ic)) &
             -(mpoly1-mpoly2)*0.5*(1.+erf((rad(k)-z2)/0.01))
      enddo
      do k=2,l
        mpoly=mpoly0-(mpoly0-mpoly1)*0.5*(1.+erf((rad(k)-z1)/width_ic)) &
             -(mpoly1-mpoly2)*0.5*(1.+erf((rad(k)-z2)/width_ic))
        hb=Rg*tempm(k-1)*((rds+zcr(k-1))/rds)**2/g0
!        tempm(k)=tempm(k-1)*(1.-(rds+zcr(k-1))*(zcr(k)-zcr(k-1))/
!     &          ((n10+1.)*hb*(rds+zcr(k))))
!        rhom(k)=rhom(k-1)*(1.-(rds+zcr(k-1))*(zcr(k)-zcr(k-1))/
!     &            ((n10+1.)*hb*(rds+zcr(k))))**(n10)
        tempm(k)=tempm(k-1)*(1.-(rds+zcr(k-1))*(zcr(k)-zcr(k-1)) &
               / ((mpoly+1.)*hb*(rds+zcr(k))))
        rhom(k)=rhom(k-1)*(1.-(rds+zcr(k-1))*(zcr(k)-zcr(k-1)) &
               / ((mpoly+1.)*hb*(rds+zcr(k))))**(mpoly)
      enddo
!
      open(1,file="strat.dat")
      do k=1,l 
        th_mlt(k)=tempm(k)*(rho00*temp00/(rhom(k)*tempm(k)))**cap
        write(1,'(F6.4,2x,F13.3,2x,F8.4,2x,F13.3)'),rad(k)/solar_radii,&
                  tempm(k),rhom(k), th_mlt(k)
        print*,rad(k)/solar_radii, tempm(k), rhom(k), th_mlt(k)
      enddo
      close(1)

    return
    end

  function deriv3(xx, xi, yi, ni, m)
!====================================================================
! Evaluate first- or second-order derivatives 
! using three-point Lagrange interpolation 
! written by: Alex Godunov (October 2009)
!--------------------------------------------------------------------
! input ...
! xx    - the abscissa at which the interpolation is to be evaluated
! xi()  - the arrays of data abscissas
! yi()  - the arrays of data ordinates
! ni - size of the arrays xi() and yi()
! m  - order of a derivative (1 or 2)
! output ...
! deriv3  - interpolated value
!============================================================================*/

implicit none
integer, parameter :: n=3
real deriv3, xx
integer ni, m
real xi(ni), yi(ni)
real x(n), f(n)
integer i, j, k, ix

! exit if too high-order derivative was needed,
if (m > 2) then
  deriv3 = 0.0
  return
end if

! if x is ouside the xi(1)-xi(ni) interval set deriv3=0.0
if (xx < xi(1) .or. xx > xi(ni)) then
  deriv3 = 0.0
  return
end if

! a binary (bisectional) search to find i so that xi(i-1) < x < xi(i)
i = 1
j = ni
do while (j > i+1)
  k = (i+j)/2
  if (xx < xi(k)) then
    j = k
  else
    i = k
  end if
end do

! shift i that will correspond to n-th order of interpolation
! the search point will be in the middle in x_i, x_i+1, x_i+2 ...
  i = i + 1 - n/2

! check boundaries: if i is ouside of the range [1, ... n] -> shift i
if (i < 1) i=1
if (i + n > ni) i=ni-n+1

!  old output to test i
!  write(*,100) xx, i
!  100 format (f10.5, I5)

! just wanted to use index i
ix = i

! initialization of f(n) and x(n)
do i=1,n
  f(i) = yi(ix+i-1)
  x(i) = xi(ix+i-1)
end do

! calculate the first-order derivative using Lagrange interpolation
if (m == 1) then
    deriv3 =          (2.0*xx - (x(2)+x(3)))*f(1)/((x(1)-x(2))*(x(1)-x(3)))
    deriv3 = deriv3 + (2.0*xx - (x(1)+x(3)))*f(2)/((x(2)-x(1))*(x(2)-x(3)))
    deriv3 = deriv3 + (2.0*xx - (x(1)+x(2)))*f(3)/((x(3)-x(1))*(x(3)-x(2)))
! calculate the second-order derivative using Lagrange interpolation
  else
    deriv3 =          2.0*f(1)/((x(1)-x(2))*(x(1)-x(3)))
    deriv3 = deriv3 + 2.0*f(2)/((x(2)-x(1))*(x(2)-x(3)))
    deriv3 = deriv3 + 2.0*f(3)/((x(3)-x(1))*(x(3)-x(2)))
end if
end function deriv3
